Method for establishing prediction model based on multidimensional texture of brain nuclear magnetic resonance images

ABSTRACT

Disclosed is a method for establishing a prediction model based on a multidimensional texture of brain nuclear magnetic resonance images. Images are segmented using a region growing method, a contourlet transform method is used to extract an edge texture feature parameter of ROIs, a multidimensional database is established, and a prediction model is established using various data mining methods, comprising a Gaussian process, a support vector machine, a random forest, a Lasso regression and a semi-supervised support vector machine. The ROIs comprise a hippocampus region and an entorhinal cortex region.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of International Patent Application No. PCT/CN2014/001129 with an international filing date of Dec. 16, 2014 designating the United States, now pending and further claims priority benefits to Chinese Patent Application No. 201410023582.5 filed Jan. 17, 2014. The content of the aforementioned applications, including any intervening amendments thereto, is incorporated herein by reference.

FIELD OF THE INVENTION

The present application relates to the field of medicine, specifically to a method for establishing a prediction model based on multidimensional textures of brain nuclear magnetic resonance images.

BACKGROUND OF THE INVENTION

In aided diagnosis of early Alzheimer disease (AD), it is important to identify features of the ROIs (reactive oxygen intermediates including entorhinal cortex and hippocampus) in the magnetic resonance imaging (MRI) images. However, the MRI technology can only use the hippocampus atrophy as one of the indicators to distinguish a patient from a healthy person. The explanations of doctors to the MRI images are easily affected by subjectivity and lack of consistency. Moreover, it is hard for the MRI technology to evaluate correctly how serious the symptom of the AD is.

1, the Prior Image Processing Technology

Contourlet Transform

The Contourlet transform inherits the anisotropic scale relationship of Curvelet transform. In a certain sense, the Contourlet transform is another implementation of the Curvelet transform. The basic idea of Contourlet transform is to capture edge point discontinuities using multiscale decomposition similar to wavelet firstly, and then to converge the point discontinuities with similar location to form contour segments according to direction information.

The Contourlet transform may be divided into two parts: Laplacian Pyramid (LP) filter structure and 2-D directional filter bank (DFB). The LP decomposition first produce a lowpass sampling approach of the original signals and a difference image of the original image and a lowpass predication image, and further decompose the lowpass image obtained to gain a lowpass image and a difference image of the next level. Multiresolution decomposition of image is obtained in this gradually filtering. The DFB adopts conjugate mirror filters with a sector structure, and transforms one level binary tree structure of a directional filter into a structure with 21 parallel passages.

The Contourlet transform is a new two-dimensional transform method for image representations, which has properties of multiresolution, localization, directionality, critical sampling and anisotropy. Its basic functions are multiscale and multidimensional. The contours of original images, which are the dominant features in natural images, can be captured effectively with a few coefficients by using Contourlet transform.

However, when the new transform method is used to process MRI images of different body parts, it needs to reconstruct new algorithm and select suitable parameters using basis functions. Therefore, there still are many theoretical issues to be studied. The Contourlet transform has been successfully used to solve actual problems such as image fusion. However, there are few references on extraction of image texture features using the Contourlet transform. No one has disclosed to use the second generation wavelet transform and Contourlet transform to extract textures of AD brain MRI images.

2, the Prior Typical Prediction Model Includes Supervised Machine Learning Model and Semi-Supervised Machine Learning Model

Supervised Machine Learning Model

Gaussian Process (GP):

The Gaussian process classification algorithm is a supervised learning algorithm, which is established under the Bayesian framework, and it is first used to solve the problem of nonlinear value prediction. Provided there is a data set D={(x_(i), y_(i))|y_(i)=±1, i=1, 2, . . . , n}, the essence of binary classification is to find a suitable mapping f(x) to ensure that the sample can be classified correctly using y=f(x). The basic idea of classification using Gaussian regression algorithm is to find a suitable f(x) aiming at maximizing the posterior probability in Bayesian rule under the assumption that the f(x) is a Gaussian process. It is firstly assumed that there is an implicit function f(x)=Ø(x)^(T)w, which defines the mapping relationship between input properties and output properties (class tags). It is further assumed that the class tags are independent identically distributed, and are subjected to conditional probability p(y|f(x)). Thus, for the binary classification, the probability that the sample x_(i) belongs to class y_(i) is: p(y _(i) |f _(i))=σ(y _(i) f _(i))

wherein, f_(i)=f(x_(i)) is an implicit function, σ(·) is a Sigmoid class function, such as a logistic function or a cumulative Gaussian function. Since training samples are independent to each other, the corresponding likelihood probability is:

${p\left( y \middle| f \right)} = {{\prod\limits_{i = 1}^{n}\;{p\left( y_{i} \middle| {f\left( x_{i} \right)} \right)}} = {\prod\limits_{i = 1}^{n}\;{\sigma\left( y_{i} \middle| f_{i} \right)}}}$

wherein, for the sample set D={(x_(i), y_(i))|y_(i)=±1, i=1, 2, . . . , n}, the Gaussian process classifier assumes that the f(x)˜GP(0, k) is a Gaussian process with mean value equivalent to zero to simplify calculation, i.e. p(f|x) is a multidimensional Gaussian density function:

${p\left( f \middle| x \right)} = {\frac{1}{\left. \left( {2\pi} \right)^{L\text{/}2} \middle| k \right|^{1\text{/}2}}\exp\left\{ {{- \frac{1}{2}}f^{T}K^{- 1}f} \right\}}$

wherein, K is a covariance matrix, which is a kernel function. Therefore, the posterior probability can be obtained according to Bayesian formula:

${p\left( {\left. f \middle| y \right.,X} \right)} = \frac{{p\left( y \middle| f \right)}{p\left( f \middle| X \right)}}{p\left( y \middle| X \right)}$

wherein, p(y|f) is a likelihood function, p(y|X)=∫p(y|f)p(f|x)df is a marginal probability density.

The training to the Gaussian process classifier is namely the estimation to the posterior probability. The estimated value q(f|y, X) of the posterior probability p(f|y,X) can be calculated by Laplace method. Therefore, the probability that the class tag of a given testing data point x* belongs to positive class is: π*=p(y*=±1|X,y,x*)=∫Ø(f*)q(f*|X,y,x*)df*

According to the principle of Gaussian process classification, the numeric data set obtained after feature extraction is inputted into the Gaussian process classifier. The estimated value q(f|y, X) of the posterior probability can be calculated after hyper-parameter optimization. Then, prediction results related to testing set can be obtained.

1) Support Vector Machine (SVM):

Support vector machine method is a machine learning method based on VC dimension theory of statistical learning theory and structure risk minimum principle. Its principle is to find an optimal classification hyperplane which satisfies the requirements of classification, so that the hyperplane ensures classification accuracy and the margin of the hyperplane is largest.

In theory, the support vector machine is capable of realizing the optimal classification of linear classifiable data. Taking the classification of two kinds of data as an example, it is provided that the training sample set is (x_(i), y_(i)), i=1, 2, . . . , l, xε{±+1}, and the hyperplane is written as (w·x_(i))+b=0. In order to ensure that the classification plane is correct to all the samples and has a classification space, the following restriction should be satisfied: y_(i)[(w·x_(i))+b]≧1, i=1, 2, . . . , l. It can be calculated that the classification space is 2/∥w∥. Thus, the problem of constructing an optimal hyperplane has been transformed to the problem of solving the following formula under the restriction formula: minØ(w)=½∥w∥ ²=½(w·w)

In order to solve this problem of restriction optimization, the Lagrange function is introduced: L(w,b,a)=½∥w∥−α(y(w·x)+b)−1

wherein, α>0, and it is a Lagrange multiplier. The solution of the problem of restriction optimization is determined by the saddle point of the Lagrange function. Moreover, the partial derivative of the solution of the problem of restriction optimization at the saddle point with respect to w and b is zero. Thus, the QP (quadratic programming) problem is transformed to a corresponding dual problem:

${\max\;{Q(\alpha)}} = {{\sum\limits_{j = 1}^{l}\;\alpha_{j}} - {\frac{1}{2}{\sum\limits_{i = 1}^{l}\;{\sum\limits_{j = 1}^{l}\;{\alpha_{i}\alpha_{j}y_{i}{y_{j}\left( {x_{i} \cdot x_{j}} \right)}}}}}}$ ${{s.t.\mspace{14mu}{\sum\limits_{j = 1}^{l}\;{\alpha_{j}y_{j}}}} = 0},{j = 1},2,\ldots\mspace{14mu},l$

the optimal solution α*=(α₁*, α₂*, . . . , α_(l)*)^(l) is obtained.

Calculating the optimal weight vector w* and optimal bias b*, which are:

$w^{*} = {\sum\limits_{j = 1}^{l}\;{\alpha_{j}^{*}y_{j}x_{j}}}$ $b^{*} = {y_{i} - {\sum\limits_{j = 1}^{l}\;{y_{j}{\alpha_{j}^{*}\left( {x_{j} \cdot x_{i}} \right)}}}}$

Wherein, the subscript jε{j|α_(j)*>0)}. Thus, the optimal classification hyperplane is obtained: (w*·x)+b*=0, and the optimal classification function is:

${{f(x)} = {{{sgn}\left\{ {\left( {w^{*} \cdot x} \right) + b^{*}} \right\}} = {{sgn}\left\{ {{\sum\limits_{j = 1}^{l}\;{\alpha_{j}^{*}{y_{j}\left( {x_{j} \cdot x_{i}} \right)}}} + b^{*}} \right\}}}},{x \in R^{n}}$

For the linear unclassifiable case, the main idea of SVM is to map the input vector to a high dimensional eigenvector space, and construct the optimal classification plane in the feature space.

Doing transform Ø to x from the input space R^(n) to the feature space H, it can be obtained: x→Ø(x)=(Ø₁(x),Ø₂(x), . . . Ø_(l)(x))^(T)

Replacing the input vector x with the eigenvector Ø(x), the optimal classification function can be obtained:

${f(x)} = {{{sgn}\left( {{w \cdot {\varnothing(x)}} + b} \right)} = {{sgn}\left( {{\sum\limits_{j = 1}^{l}\;{\alpha_{i}y_{i}{\varnothing(x)}}} + b} \right)}}$

2) Random Forests (RF):

Random forest is a classification and predication model proposed by Leo Breiman, which has been applied in many fields.

The basic idea of constructing random forests is continuously generating training sample and testing sample by the technique of bootstrap resampling. The training sample generates a plurality of classification tree constituting the random forests. The classification result of the testing data is determined according to the vote of classification tree.

It is provided that the content of original sample is N, the variables representing various genes are x₁, x₂, . . . , x_(m). Randomly sampling with replacement new self-service samples with the number of b, and hence forming classification trees with the number of b. The samples that are not sampled each time form the data out of bag (OOB) with the number of b. the data out of bag acting as testing sample can be used to evaluate the importance of the various variables in the classification. The realization process is as follows:

Using the self-service samples to form each tree classifier, and classifying the corresponding OOB at the same time, it can hence obtain the voting scores of each sample in the OOB of the self-service samples with the number of b, which can be written as vote₁, vote₂, . . . , votes.

Randomly changing the order of the value of the variables x_(i) in the OOB samples with the number of b, it forms new OOB testing samples. Then, classifying the new OOB by the constructed random forests and obtaining the voting scores of each sample according to the number of sample that is judged to be correct.

3) Lasso:

The basic idea of Lasso (the least absolute shrinkage and selection operator) method is to minimize the sum of squared residuals under the restriction that the sum of absolute regression coefficient is less than a constant, thus generating some regression coefficient rigidly equaling to zero. Although this is a biased estimation to decrease the variance of the predicted value, the prediction accuracy of the model will be increased accordingly and the model obtained is easier to be explained. For the case with a plurality of independent variables, we always hope to determine a model with fewer variables to present the best effect. The selection criteria of Lasso is to minimize the penalty function

${{\frac{1}{n}{\sum\limits_{i = 1}^{n}\;\left( {Y_{i} - {X_{i}\beta_{i}}} \right)^{2}}} + {\sum\limits_{j = 1}^{\lambda}\;{\beta_{j}}}},$ and adopts the least-angle regression algorithm.

The regression steps of Lasso are as follows:

a) Provided the estimated initial value of β is β⁽⁰⁾, fixing λ and selecting the least-squares estimation {tilde over (β)}_(LS) as the initial value of β. That is to say, β⁽⁰⁾={tilde over (β)}_(LS)=(X′X)⁻¹X′Y;

b) Obtaining the estimation to β(written as {tilde over (β)}) when fixing λ according to the estimated initial value β⁽⁰⁾ and the following iterative formula: β⁽¹⁾)=β⁽⁰⁾−[∇² Q _(LS)β⁽⁰⁾ +nΣ _(λ)β⁽⁰⁾]⁻¹[∇Q _(LS)β⁽⁰⁾ +nU _(λ)β⁽⁰⁾]

wherein, ∇²Q_(LS)β⁽⁰⁾=X′X, ∇Q_(LS)β⁽⁰⁾=−2X(Y−Xβ⁽⁰⁾), Σ_(λ)β⁽⁰⁾=λdiag(1/|β₀|, . . . 1/|β_(p)|), if β_(j)=0, the value of 1/|β_(j)| can be evaluated to be a very large number. U_(λ)β⁽⁰⁾=U_(λ)β⁽⁰⁾β⁽⁰⁾;

-   -   c) Changing the value of λ in the range of 0˜1, under the         restriction of λ, constructing generalized cross validation

${{GCV}(\lambda)} = {\frac{{{Y - {X\;{\overset{\sim}{\beta}}_{\lambda}}}}^{2}}{\left( {1 - {{e(\lambda)}/n}} \right)^{2}}.}$ The penalty function that minimizes the cross validation is the optimal selection penalty parameter λ, λ=arg min GCV(λ), and hence the Lasso estimation to β is {circumflex over (β)}_(Lasso)={circumflex over (β)}_(Lasso).

-   -   (1) Semi-Supervised Support Vector Machine (S3VM):

The semi-supervised learning is a key problem studied in the field of pattern recognition and machine learning in recent years, which mainly considers the problem of how to train and classify with a few tag samples and a plurality of samples without tag. The semi-supervised learning algorithm based on support vector machine is a model first proposed by Kristin, which is realized by mixed integer programming method. Later, Joachims, Demirez et al. proposed some similar methods one after another among which the transductive support vector machine (TSVM) algorithm is more typical.

Given that a training set T={(x₁, y₁), . . . , (x_(l), y_(l)), x_(j+1), . . . , x_(j+k)}, wherein, x_(f)εX=R_(n), the top l samples are samples with classification tag, i.e. when i=1, 2, . . . , l, y_(i)εY={−1,1}; the down k samples are samples without classification tag. Finding a decision function f(x)=sgn(ω^(T)φ(x)+b) at X=R^(n) to deduce the class (positive class or negative class) of x of any pattern (wherein, ω is weight vector, φ(·) is a mapping function, b is a constant, and the kernel function K(x_(i), x_(j))=<φ(x_(i)), φ(x_(j))>). Accordingly, the solving process is actually to find a principle to divide the points in R^(n) into two parts.

Obviously, it exists ωεR^(n), bεR, to ensure that it exists y_(i)(ω^(T)φ(x_(i))+b)+τ_(i)≧1, ξ_(i)≧0, i=1, . . . , l for each sample with classification tag x_(i)(i=1, . . . , l). It exists (ω^(T)φ(x_(j))+b)+r_(j)≧1, ω^(T)φ(x_(j))+b−s_(j)≦−1, r_(j), s_(j)≧0, j=l+1, . . . , l+k.

Currently, such a classification learning machine is constructed by maximizing soft margin and slack variables ξ=(ξ₁, ξ₂, . . . , ξ_(l)), r=(r_(l+1), r_(l+2), . . . , r_(l+k)), s=(s_(l+1), s_(l+1), . . . , s_(l+k)) are introduced. Obviously, ξ, r, s reflect the case that the training set is wrongly classified. Thus, it can construct the degree that describes the training set wrongly classified, and it may as well to adopt Σ_(i=1) ^(l)ξ_(i) ², Σ_(j=1) ^(l+k)r_(j)s_(i) as a mensuration to degree that describes the training set wrongly classified. At the same time, it is expected that the wrongly classified degree Σ_(i=1) ^(l)ξ_(i) ², Σ_(j=1) ^(l+k)r_(j)s_(i) is as small as better, that the optimal classified hyperplane can be solved out by minimizing ½(ω^(T)ω+b²), and that the affection to the classification result is limited. In order to fully take the three aims into account, the penalty parameters c₁ and c₂ are introduced as a weight to balance the three aims. That is to say, the aim function is minimized as:

${\frac{1}{2}\left( {{\omega^{T}\omega} + b^{2}} \right)} + {c_{1}{\sum\limits_{i = 1}^{l}\;\xi_{i}^{2}}} + {c_{2}{\sum\limits_{j = 1}^{l + k}\;{r_{j}s_{i}}}}$

Thus, the problem of solving the optimal classified hyperplane is transformed into the problem of optimization:

$\mspace{20mu}{\min\limits_{\omega,b,\xi,r,s}\left( {{\frac{1}{2}\left( {{\omega^{T}\omega} + b^{2}} \right)} + {c_{1}{\sum\limits_{i = 1}^{l}\;\xi_{i}^{2}}} + {c_{2}{\sum\limits_{j = 1}^{l + k}\;{r_{j}s_{i}}}}} \right)}$   s.t.  y_(i)(ω^(T)φ(x_(i)) + b) + ξ_(i) ≥ 1, i = 1, …  , l   ω^(T)φ(x_(j)) + b + r_(j) ≥ 1, ω^(T)φ(x_(j)) + b − s_(j) ≤ −1, r_(j), s_(j) ≥ 0,  j = l + 1, …  , l + k   Given ${L\left( {\omega,b,\xi,r,\alpha,\beta,\gamma,u,v} \right)} = {{\frac{1}{2}\left( {{\omega^{T}\omega} + b^{2}} \right)} + {c_{1}{\sum\limits_{i = 1}^{l}\;\xi_{i}^{2}}} + {c_{2}{\sum\limits_{j = 1}^{l + k}\;{r_{j}s_{i}}}} - {\sum\limits_{i = 1}^{l}\;{\alpha_{i}\left\lbrack {{y_{i}\left( {{\omega^{T}{\varphi\left( x_{i} \right)}} + b} \right)} + \xi_{i} - 1} \right\rbrack}} - {\sum\limits_{j = {l + 1}}^{l + k}\;{\beta_{j}\left( {{\omega^{T}{\varphi\left( x_{j} \right)}} + b + r_{j} - 1} \right)}} - {\sum\limits_{j = {l + 1}}^{l + k}\;{\gamma_{j}\left( {{{- \omega^{T}}{\varphi\left( x_{j} \right)}} - b + r_{j} - 1} \right)}} - {\sum\limits_{j = {l + 1}}^{l + k}\;{u_{j}r_{j}}} - {\sum\limits_{j = {l + 1}}^{l + k}\;{v_{j}s_{j}}}}$ $\mspace{20mu}{{Thus},{\frac{\partial L}{\partial\omega} = {0 = {{> \omega} = {{\sum\limits_{i = 1}^{l}\;{y_{i}\alpha_{i}{\varphi\left( x_{i} \right)}}} + {\sum\limits_{j = {l + 1}}^{l + k}\;{\left( {\beta_{j} - \gamma_{j}} \right){\varphi\left( x_{j} \right)}}}}}}},\mspace{20mu}{\frac{\partial L}{\partial b} = {0 = {{> b} = {{\sum\limits_{i = 1}^{l}\;{y_{i}\alpha_{i}}} + {\sum\limits_{j = {l + 1}}^{l + k}\;\beta_{j}} - \gamma_{j}}}}},\mspace{20mu}{\frac{\partial L}{\partial b} = {0 = {{> \xi_{i}} = {\frac{1}{2\; c_{1}}\alpha_{i}}}}},\mspace{20mu}{\frac{\partial L}{\partial b} = {0 = {{> s_{j}} = \frac{\beta_{j} + u_{j}}{c_{2}}}}},\mspace{20mu}{\frac{\partial L}{\partial b} = {0 = {{> r_{j}} = \frac{\gamma_{j} + v_{j}}{c_{2}}}}},}$

It can be obtained with substitution:

$\min\limits_{\omega,b,\xi,r,s}\left( {{\frac{1}{2}{\sum\limits_{i_{1} = 1}^{l}\;{\sum\limits_{i_{2} = 1}^{l}\;{y_{i_{1}}y_{i_{2}}\alpha_{i_{1}}\alpha_{i_{2}}{K\left( {x_{i_{1}},x_{i_{2}}} \right)}}}}} + {\frac{1}{2\; c_{1}}{\sum\limits_{i = 1}^{l}\;\alpha_{i}^{2}}} + {\frac{1}{2}\left( {{\sum\limits_{i = 1}^{l}\;{y_{i}\alpha_{i}}} + {\sum\limits_{j = {l + 1}}^{l + k}\;\left( {\beta_{j} - \gamma_{j}} \right)}} \right)^{2}} + {\frac{1}{2}{\sum\limits_{j_{1} = {l + 1}}^{l + k}\;{\sum\limits_{j_{2} = {l + 1}}^{l + k}\;{\left( {\beta_{j_{1}} - \gamma_{j_{1}}} \right)\left( {\beta_{j_{2}} - \gamma_{j_{2}}} \right){K\left( {x_{j_{1}},x_{j_{2}}} \right)}}}}} + {\sum\limits_{i = 1}^{l}\;{\sum\limits_{j = {l + 1}}^{l + k}\;{y_{i}{\alpha_{i}\left( {\beta_{j} - \gamma_{i}} \right)}{K\left( {x_{i},x_{j}} \right)}}}} + {\frac{1}{c_{2}}{\sum\limits_{j = {1 + k}}^{l + k}\;{\left( {\beta_{j} + u_{j}} \right)\left( {\gamma_{j} + v_{j}} \right)}}} - {\sum\limits_{i = 1}^{l}\;\alpha_{i}} - {\sum\limits_{j = {1 + k}}^{l + k}\;\left( {\beta_{j} + \gamma_{j}} \right)}} \right)$   s.t.  α_(i), β_(j), γ_(j), u_(j), v_(j) ≥ 0.

SUMMARY OF THE INVENTION

The objective of the present invention is to provide a method for establishing prediction models based on multi-dimensional textures of brain nuclear magnetic resonance images so as to establishing various kinds of prediction models in order to effectively find the MCI patient, diagnose early AD and observe the brain structure change of healthy geriatric population.

The present invention synthesizes the advantages of region growing method in the prior art, utilizes Contourlet transform to extract margin texture features of ROIs. The texture extraction method is comprehensive and novel, which ensures the realization of the objective of the present invention.

The present invention adopts the following technical route to establish a method for segmenting MRI images of ROIs of geriatric population who has MCI and early AD nidus and who are healthy, and for extracting texture features (specifically see FIG. 3), and prediction models for judging related ROIs features have been established accordingly:

1) establishing a ROIs database;

2) segmenting the ROIs in the images using the region growing method;

3) processing the images using the Contourlet transform to extract the following variables: entropy, mean gray, correlation, energy, homogeneity, variance, maximum probability, inverse difference moment, cluster tendency, contrast, sum mean, difference mean, sum entropy and difference entropy;

4) establishing a database of image feature parameters using the variables obtained in steps (2) and (3);

5) establishing the prediction model according to the database in step (4), wherein methods for establishing the prediction model comprises Gaussian process, support vector machine, random forests and Lasso regression; and

6) repeatedly verifying the feature parameters obtained in step (4) using a sample to modify the prediction model so as to obtain an ideal model with correct results.

At the present, there is no uniform definition to the textures of images. It is widely accepted that the textures of images describe the surface gray of objects or the change of color. Such change is related to the properties of the object, which is a repetition of a certain texture primitive.

Adopting Contourlet transform, the following texture feature parameters can be obtained:

Entropy: f₄=Σ_(i=0) ^(K−1)Σ_(j=0) ^(K−1)c_(ij) log c_(ij) reflects the randomness of Image's texture;

Mean gray:

$f_{15} = {\frac{1}{Q}\Sigma\;{p\left( {i,j} \right)}}$ reflects the gathering tendency of all the gray values of the image;

Correlation:

$f_{5} = \frac{{\sum\limits_{i = 0}^{K - 1}\;{\sum\limits_{j = 0}^{K - 1}\;{({ij})c_{ij}}}} - {u_{x}u_{y}}}{\sigma_{x}\sigma_{y}}$ measures the correlation of the image gray;

Energy (angular second moment): f₁=Σ_(i=0) ^(K−1)Σ_(j=0) ^(K−1)c_(ij) ² reflects distribution uniformity of image gray and thickness of the texture;

Homogeneity:

$f_{13} = {\sum\limits_{i = 1}^{M}\;{\sum\limits_{j = 1}^{N}\;\left( \frac{P\left( {i,j} \right)}{1 + {{i - j}}} \right)}}$ reflects the similarity of the gray;

Variance:

$f_{11} = {\frac{1}{Q - 1}{\sum\limits_{i = 0}^{N - 1}\;\left\lbrack {{P\left( {i,j} \right)} - \mu} \right\rbrack^{2}}}$ reflects the distribution of the gray;

Maximum probability (MP):

$f_{514} = \frac{{frequency}_{\max}}{Q}$ is the incidence of most outstanding pixel pairs;

Inverse difference moment (IDM):

$f_{3} = {\sum\limits_{i = 0}^{K - 1}\;{\sum\limits_{j = 0}^{K - 1}\;{\frac{1}{1 + \left( {i + j} \right)^{2}}C_{ij}}}}$ reflects smoothness of the image;

Cluster tendency: f₁₀=Σ_(i) ^(M)Σ_(j) ^(N)(i−j)² P(i,j) reflects the resolution of the image;

Sum mean: f₆=Σ_(k=0) ^(2k−2)kc_(x+y)(k) and difference mean: f₇=Σ_(k=0) ^(k−1)kc_(x−y)(k) provide the mean value of the image gray;

Sum entropy: f₈=Σ_(k=0) ^(2k−2)c_(x+y)(k)log {c_(x+y)(k)} and difference entropy: f₉=Σ_(k=0) ^(k−1)c_(x−y)(k)log {c_(x−y)(k)}.

Adopting above method of the present invention to establish prediction models for early AD, the judging accuracy is increased to 100%.

An example of extracting textures in the ROIs in provided below, its steps are as follows:

1, collecting original MRI images of the elderly of AD, MCI and the healthy elderly (the MRI image of the healthy elderly is provided in the drawings for example), as shown in FIG. 1;

2, segmenting the above images by region growing method to obtain images as shown in FIGS. 2 and 3; the segmentation is executed by directly running a written program;

3, extracting the texture feature parameters utilizing Contourlet transform, there is a corresponding program for each Contourlet transform method and the Contourlet transform can be executed by directly running the program.

Advantageous Effects

It is one innovation of the present invention to apply the region growing method to segment brain MRI images. It is proved by experiment that extracting texture by the region growing method to model is better than extracting texture by overall segmenting to model, which can reserve the margin information better.

Contourlet transform is an effective method to extract texture features in the brain MRI images, which can extract 14 texture feature parameters so as to reflect the texture features of images more fully.

It is proved by experiment that the specificity for predicting the sensitivity of early AD is 100% by extracting texture feature parameters from ROIs of brain MRI images of small sample to establish support vector machine prediction model. The small sample is obtained by utilizing the region growing method to segment ROIs, selecting 80% data as training samples and the rest 20% for verification samples, so as to test the consistency between determination of ROIs pathology by utilizing model and determination of ROIs by adopting pathological diagnosis.

From above data, it can be concluded that the support vector machine prediction model, which is established by adopting region growing method to segment ROIs of the brain MRI images and extracting texture feature parameters by Contourlet transform, takes a good effect in early AD added diagnosis.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is an original MRI image of a healthy elderly, wherein ROIs are marked with black boxes;

FIG. 2 is an enlarged image of left brain ROIs which is obtained by segmenting FIG. 1 by region growing method;

FIG. 3 is an enlarged image of right brain ROIs which is obtained by segmenting FIG. 1 by region growing method;

FIG. 4 is technical route for extracting multi-dimensional textures of the brain MRI images using the method of the present invention;

In the accompanying drawings: 1, left brain hippocampus, entorhinal cortex region; 2, right brain hippocampus, entorhinal cortex region.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

The following embodiments are used to describe the model established by using the method of the present invention for predicting AD based on brain MRI images including related ROIs. They are only used to further explain the method of the present invention. However, the embodiments are not used to limit the application scope of the present invention. Actually, the method can also be used to determine properties of other type medical images.

The source of the images: brain MRI images of the elderly of AD, MCI and the healthy elderly shared on ADNI website, which are in form of .Nii and are read by MRIcro software.

Methods: programming with Matlab software, segmenting the ROIs in above MRI images with region growing method, and extracting texture feature parameters of related ROIs with Contourlet transform.

An embodiment for extracting texture feature parameters ROIs of brain MRI images is provided below; its steps are as follows:

1, respectively collecting a total of 504 original brain MRI images from 15 cases of AD illness, 15 cases of MCI illness, and 15 cases of healthy persons.

2, segmenting the images by the region growing method. A general method of the region growing method is used to segment the images. A programmed procedure is adopted, the threshold is set to be 35, and the procedure is directly executed. FIG. 1 is an example of image to be segmented and FIGS. 2 and 3 are figures obtained after segmentation.

3, adopting Contourlet transform to respectively extract the texture feature parameters.

4, adopting Contourlet transform to extract ROIs of brain MRI images to obtain 14 texture feature parameters to establish a prediction model.

When using the support vector machine model to predict the AD group and the MCI group, the classification effects are that the sensitivity is 100% and the specificity is 100% respectively. When predicting the AD group and the healthy group, the classification effects are that the sensitivity is 100% and the specificity is 100% respectively. When using the Lasso model to predict the AD group and the MCI group, the classification effects are that the sensitivity is 83.3% and the specificity is 100% respectively. When predicting the AD group and the healthy group, the classification effects are that the sensitivity is 73.7% and the specificity is 87.5% respectively.

It is proved by experiment that the specificity for predicting the sensitivity of early AD is 100% by extracting texture feature parameters from ROIs of brain MRI images of small sample to establish support vector machine prediction model. The small sample is obtained by utilizing the region growing method to segment ROIs, selecting 80% data as training samples and the rest 20% for verification samples, so as to test the consistency between determination of ROIs pathology by utilizing model and determination of ROIs by adopting pathological diagnosis.

From above data, it can be concluded that the support vector machine prediction model, which is established by adopting region growing method to segment ROIs of the brain MRI images and extracting texture feature parameters by Contourlet transform, takes a good effect in early AD added diagnosis. 

We claim:
 1. A method for establishing a prediction model based on multidimensional textures of brain nuclear magnetic resonance images, comprising segmenting the images using a region growing method, extracting, edge texture feature parameters of ROIs using Contourlet transform, establishing a database, and establishing the prediction model, the method specifically comprising the following steps: (1) establishing a ROIs database; (2) segmenting the ROIs in the images using the region growing method; (3) processing the images using the Contourlet transform to extract the following variables: entropy, mean gray, correlation, energy, homogeneity, variance, maximum probability, inverse difference moment, cluster tendency, contrast, sum mean, difference mean, sum entropy and difference entropy; (4) establishing a database of image feature parameters using the variables obtained in steps (2) and (3); (5) establishing the prediction model according to the database in step (4), wherein methods for establishing the prediction model comprises Gaussian process, support vector machine, random forests, Lasso regression and semi-supervised support vector machine; and (6) verifying the feature parameters obtained in step (4) using a sample to modify the prediction model so as to obtain an ideal model with correct results. 